Unified continuum approach to crystal surface morphological relaxation 
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A continuum theory is used to predict scaling laws for the morphological relaxation of crystal 
surfaces in two independent space dimensions. The goal is to unify previously disconnected 
experimental observations of decaying surface profiles. The continuum description is derived from 
the motion of interacting atomic steps. For isotropic diffusion of adatoms across each terrace, 
induced adatom fluxes transverse and parallel to step edges obey different laws, yielding a tensor 
mobility for the continuum surface flux. The partial difi'erential equation (PDF) for the height 
profile expresses an interplay of step energetics and kinetics, and aspect ratio of surface topography 
that plausibly unifies observations of decaying bidirectional surface corrugations. The PDF reduces 
to known evolution equations for axisymmetric mounds and one-dimensional periodic corrugations. 



Novel small devices rely on the stability of nanoscale 
surface features. The lifetimes of nanostructures decay- 
ing via surface diffusion scale as a large power of their size 
and increase with decreasing temperature. Below rough- 
ening, crystal surfaces evolve via the motion of atomic 
steps bounding nanoscale terraces [TJ [5] . 

Experiments with decaying surface features O SI [5l |6j 
III HI H] are useful for testing step models. Particularly 
informative are observations of bidirectional corrugations 
relaxing below roughening [51 [SJ |7J HI IS] ■ In lithography- 
based experiments [5], where initial wavelengths in two 
directions differ significantly and profiles depend nearly 
on one space dimension (ID), the surface height decays 
exponentially with time. By contrast, in sputter-rippling 
experiments [71 |S], where initial wavelength ratios are 
closer to unity and profiles evidently depend on two space 
dimensions (2D), height spatial- frequency components 
decay inverse linearly with time. These observations have 
previously evaded a unified theory [TUl |TT] . In this Let- 
ter, we use a continuum theory to plausibly unify these 
observations via an appropriate tensor mobility. 

There are two main theoretical approaches to crystal 
surface morphological evolution below roughening. One 
approach follows the motion of steps in the spirit of the 
Burton-Cabrera-Frank (BCF) model [T] via numerical so- 
lutions of coupled equations for step positions [12 IT^ . 
Step simulations in ID [12] show exponential decay of 
surface corrugations with attachment-detachment lim- 
ited (ADL) kinetics, in agreement with lithography ex- 
periments [5]. Step simulations in 2D invoke axisymme- 
try |13j . and are thus limited in their ability to make 
predictions for general surface morphologies. 

Another approach relies on equilibrium thermodynam- 
ics and mass conservation using continuum evolution 
laws [IldllllSlinillHlIinillD] such as partial differential 
equations (PDEs). PDEs enable simple scaling predic- 
tions; see e.g. [17]. Continuum models are criticized [12] 
for their inaccurate description of macroscopic, planar 



surface regions ( "facets" ) , but progress is made in includ- 
ing facets in evolution laws |21j . Continuum theories have 
not previously unified observations of decaying surface 
corrugations [TT]. An ingredient of such theories is the 
scalar mobihty for the adatom flux in 2D [rUllTTl [T^ . 
which does not essentially distinguish adatom fluxes par- 
allel to steps from fluxes transverse to steps. This for- 
mulation is valid when steps are everywhere parallel |17j , 
but is shown here to be inadequate in general cases. 

In this Letter we plausibly unify experimental obser- 
vations of decaying proflles by invoking a tensor macro- 
scopic mobility for the adatom flux in a setting of 
isotropic terrace diffusion; see Eqs. ([8|-(10l. An elabo- 
rate derivation is given elsewhere |22j . Here, we provide a 
more general yet simpler derivation. We show that the re- 
sulting PDE for the height proflle reduces to known evo- 
lution laws for ID gratings and 2D nanostructures. Fur- 
ther, we relate scaling predictions of the general theory 
to relaxation experiments. We find that observed scaling 
laws with time can arise from competition of step kinetics 
with surface topography. This effect is due to coupling of 
adatom flux components via terrace diffusion, and is dis- 
tinct from the influence of step edge diffusion, e.g. work 
in [23j . A similar effect of anisotropic terrace diffusion on 
step meandering is studied in [53]. By contrast to |24j . 
our step model has scalar microscopic parameters. 

First, we describe the model of step flow [1]. A top 
terrace is surrounded by non-self-intersecting and non- 
crossing steps numbered i = 1,2,...; i — 1 denotes 
the top step. The projection of steps on the basal 
(high-symmetry) plane is described by the position vec- 
tor v{ri,a,t) where t is time; r/ = r/i at the ith step, 
rji < T] < 77i_|_i on the zth terrace, and a gives the position 
along each step; see Fig. 1. The unit vectors normal and 
parallel to steps in the direction of increasing 77 and a 
are e,, and e^-; • = 0. The metric coefficients (to be 
used below) are = |9^r| and — l^^rj; dj^ := d/d-q. 

Mass conservation for atoms is described by 
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FIG. 1: Schematic of steps on the basal plane. Local coor- 
dinates relative to a top terrace are {rj,a). The arrow shows 
longitudinal flux directed to a valley. Dots denote many steps. 



Vi is the (normal) velocity of the zth step, f2 is the 
atomic volume, and a is the step height; J-' ~ 3i ■ e,, 
is the adatom current (adatoms/length/time) transverse 
to steps; Ji — —Ds^Ci is the adatom current on the zth 
terrace, Ds is the terrace difFusivity, a scalar function of 
r, and Ci(r,i) is the adatom density [adatoms/ (length)'^] 
on the ith terrace. The variable Ci solves the diffu- 
sion equation, which in the quasistatic approximation 
becomes V^C^ « 0, where no material is deposited from 
above. The requisite boundary conditions describe atom 
attachment-detachment at the ith and (z+l)th steps [T^ . 



TJ'^{vu<y')^k[C,{m,<y')-Cl\a')]. 



(2) 



The time {t-) dependence is omitted, I = i (upper sign) 
or i -|- 1 (lower sign), k is the attachment-detachment 
rate, and C'p^ is the ith step equilibrium atom density. 
Note that Eqs. ^ are similar to those appearing in 
other growth problems; but in the present case there is 
no morphological instability. 

Next, we close Eqs. ([T|) and ^ by relating C'^'^ with 
the step positions. First, we introduce the step chemical 
potential of the ith step, fj,i{a,t), the change in the step 
energy by adding or removing an atom at (rji^a) |13j : 

Cs is the atom equilibrium density near a straight isolated 
step and fceT is the Boltzmann energy. 

Second, we provide a relation of with the step posi- 
tions. We use U{ri,a), the energy of atoms per length of 
the ith step (for t] = rji); thus, the length 5si — S,a5a of 
the ith step has energy SWi = U{rii,a)Ssi. Addition or 
removal of atoms at {r]i, a) causes rji to change by Sr] as- 
suming energy isotropy, the step to move along the local 
normal (e^) by distance Sg — ^.r/Sri, and the step energy 
to change by S^Wi — [d,j.{SWi)]Sr]. By definition of fii, 
? 07 for (fo,577)->0, wefind 



M. = (r!/a)P,C/)e^-i + «C/(r/,(T) 



(3) 



where k is the step curvature and U — "f + ; 7 is the 



step line tension, assumed a constant, and C/™* accounts 
for interactions with other steps. For nearest-neighbor 
elastic-dipole or entropic repulsions, C/'"* 



U 



int 



$(r7,r7i+i;cr) <^{r],r]^ 



IS [21 [25] 
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where g (energy/length) is positive, and ^{r], (; a) is 
geometry-dependent, differentiable with (77, (^) and sat- 
isfies ^{r]i,r]i+i)Ssi = ^{r]^+i,r]i)5si+i [22]. Suppressing 
77i_i and r]i+i, Eqs. ^ and Q yield fii = jl{ri = 77^, cr, t). 

Equations ([l])-(|4]) describe coupled step motion via 
adatom isotropic diffusion across terraces and atom 
attachment-detachment at steps. To enable predictions 
for decaying surface profiles at length scales large com- 
pared to the terrace width, Sgi, we next derive a PDE for 
the continuum height profile, h(r,t). Thus, Sgi is small 
compared to: (i) the length over which the step density, 
j^, varies; and (ii) the step radius of curvature, 1/k. We 
take Srji = fji+i ~ ~* ^ with fixed In this limit, 
|V/i|, where Vh = {dxh, dyh), and Vi — > 
First, we note that the familiar continuum mass con- 
servation statement for atoms comes from the step veloc- 
ity law, Eq. ([T|). By using the continuum surface current 
3(r,t), the continuous extension of 3i{rii,a,t), we obtain 



5p 



dth- 



-[9„(c. J") + 9.(6, J")] = -^^v • J. (5) 



Next, we apply Eqs. Q to relate J(r,i) to the contin- 
uum step chemical potential, /i(r,t) = il(rji,a,t). The 
following procedure is more general than the analysis 
in [22 . (i) We apply Eq. ([2| with the upper sign for cr' — 
a, and with the lower sign for a' = a+5a. (ii) We expand 
the transverse current, J]' , the density Ci and jl, each 
evaluated at {rji+i, a+6a), at {rji, a) using — —DgVCi, 

e.g., au+i,^+sa «^ - D-\i,,5n,f^\, + iMJt\i) 

where Q\p := Q{r]p, a) and = -eg. is the longitudinal 
current, (iii) We subtract Eqs. Q dropping terms that 
are negligible as Srji — > 0. Thus, we find 



1 + -Z— 
where q = 



IDs. 
ka 



6, 



5CT 



J(r,t) 



By setting 5a = 



in Eq. ^ we obtain 
1 djjp. 



0, 
(6) 



ksT l + q\Vh\ 



where q{a/Spi) is fixed. Hence, Eq. ^ reduces to 

knT 



(7a) 



(7b) 



By Eq. (7b I, the continuum longitudinal current, J" 



has the terrace diffusivity, Dg, whereas the normal cur- 



rent J'', Eq. (7a I, has the slope-dependent effective dif- 
fusivity Ds — Ds{l + q\\/h\)~^] Ds equals Ds for terrace- 
diffusion limited (TDL) kinetics, g|V/i| ^ 1. This behav- 
ior results from coarse-graining in 2D, combining atom 



3 



attachment-detachment, terrace diffusion and step to- 
pography. For ADL kinetics, g|V/i| » 1, is sensitive 
to step variations of n because steps are sources and sinks 
of atoms by Eqs. (|2]), whereas J'^ is sensitive to space vari- 
ations of /i along steps due to adatom diffusion between 
non-parallel steps. Equations ([t]) read J = — CgM • 
where the mobility M (length^ /energy/time) is a second- 
rank tensor (a 2 x 2 matrix where J and V/i are 2-column 
vectors). In the basal's plane Cartesian system (x, y) the 
matrix elements Af^ (i, j = x, y) of M are 



kBT |V/i|2 



M, 



hq\Vh\ 
q\\/h\ (d^h) 



ksT 1 
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Myy = 



kBT |V/i|2 



l + q\Vh 



q\Vh\ |V/l|2 

,2 

1 



a. 



(8) 
(9) 
(10) 



where a :— 

X 



g-j^. For biperiodic profiles, a is estimated by 
^ , the (aspect) ratio of dominant (maximum-amplitude) 
wavelengths in x and y; we take < Xy and, thus, a < 1. 

Next, we obtain a PDE for the height profile, h{r,t). 
First, we derive a relation of fi with Vh via Eqs. ([s]) and 
Q. (i) We expand in {rji — () the function ^{rfi, (; a) of 
Eq. Q, where C = 77^+1 or rji-i. (ii) We use an identity 
for (f>, which stems from the definition of t/™' [22 . After 
some algebra, the limit 6rii yields 



(11) 



where k 



V • T^rn- is the step edge curvature, gi = ^ 



\Vh\ 

and (73 — ^{^)^^{r]i,rii); gi and g^ have dimensions 
energy per area. This /i also results from the variational 
derivative of the surface energy E = J J dxdy [gi\\7h\ + 
(53/3)|V/i|3] [TZl^. By Eqs. (§, (0 and 



dfh^ BV-i A-V 



93 
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V-(|V/i|V/i) 



(12) 



where A = and B = [(length) Vtime]. 

By Eqs. ([8])-([l0j) for M, Eq. Q describes an interplay 
of step energetics and kinetics, and aspect ratio a. This 
dependence on a is absent in previous studies of mor- 
phological evolution below roughening (TUl |T71 [HI . 

It is tempting to compare Eq. ( |12[ ) and its ingredients 
to similar continuum laws for steps, e.g. Eq. (14) of |[24j 
for a step meander without deposition. The last term of 
Eq. (14) in [21] pertains to the flux along the step edge, 
with a mobility that depends on the step edge slope. In 
the small slope limit, this term appears to agree with 
Eq. (7b I. We emphasize that the isotropic physics of 



our model is different from that of |24j where anisotropic 
terrace diffusion coexists with step edge diffusion. 

We now show that Eq. ( 12 ) reduces properly to 



First, we have J'^ = by which the effective mobility 
becomes M = ^^^(1 + 9|Vft.|)^^, a scalar. For straight 
steps (in ID), r] = x, we have k = and the PDE 
becomes dth = -B^d^Kl + q\d^h\)-^d^^{\d^h\d^h)] 
where i?3 = '"'fi^^ , which is consistent, for exam- 
pie, with [12]. The reduced PDE can be applied to 
systems of periodic corrugations in ID [S] [T5] . For 
concentric circular, descending steps in 2D, 77 cx r (polar 
distance), we have k = 1/r and the PDE ( |T2| becomes 
dth = Br-^ driil + qm)-^[-r-^ + i^rdrir^r{rm^))]} 
where m — \drh\, which is applied to decaying axisym- 
metric mounds [t] [T3l ITT] [2T] . 



We now apply separation of variables to Eq. ( 12 1 for 



smooth regions, aiming to unify decay laws in relaxation 
experiments. Consistent with step simulations in ID [12] 
and kinetic Monte Carlo simulations in 2D [19], both 
for initial sinusoidal profiles, we set h{r,t) w A{t)H{r) 
and find A(t). This variable separation, which we call 
a "scaling Ansatz" , is satisfied only approximately: 
additive terms in /i and M scale differently with A. In 
/i, Eq. (Ill, the step line tension (171 term) scales with 



A" and the step interaction (53 term) scales with A^; in 



M,j, Eqs 



-([10|, the kinetic term l3 = {l + q{\Wh\))-^ 
must be compared to the aspect ratio squared, a^; 
(|V/i|) = J/ is a typical slope. 

Our analysis does not address the evaluation of H{r), 
which solves a nonlinear PDE. Because boundary con- 
ditions for H at facet edges require feedback from step 
simulations [51], a viable numerical scheme for H is not 
possible at the moment. By [T^ldS], the scaling Ansatz 
seems reasonable for long t and initial sinusoidal profiles. 

We next focus on ADL kinetics, (3 <^ 1, distinguishing 
four cases. In the first case: (i) step interactions domi- 
nate, I53V • (|V/i|V/i)| » M or | » (^)2 by dimen- 

sional analysis for sinusoidal profiles, where v « and 
hpv is the peak-to- valley height variation; and (ii) /? < 
so that longitudinal fiuxes are considerable. Thus, fi 
scales with A'^, and the matrix elements of A are A^x ~ 

r>y-tA A A ~ 



which scale with A° as in TDL kinetics. We find A « 
—cB^A^, where the dot denotes time derivative. Hence, 

A{t) = Ao (l + cBsAot)-'; Ao := A{0). (13) 

The constant parameter c [(length)^"*] depends on H{r) 
and is thus affected by facet evolution. Equation (13 1 



A„ 



suggests that surface relaxation is inverse linear with 
time if the (y-) adatom flux in the direction of the longer 
wavelength (Ay) is significant. 

In the second case: (i) step interactions remain dom- 
inant, and (ii) P > a^, so that transverse fluxes prevail. 
Thus, we obtain A = —CB^A, by which 



A{t) = Ao e 



-CBst 



(14) 



known macroscopic laws for everywhere parallel steps. 



where C is affected by H. The remaining cases for ADL 
kinetics follow similarly. The results are summarized 
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in Table I. The square-root decay with time when line 
tension dominates and /? > is in agreement with |19j . 



TABLE I: Summary of decay laws for the height amplitude 
A{t) in ADL kinetics, j3 <^ 1. The top row lists the possible 
kinetic-geometric conditions on the mobility, Eqs. (|8|-(10l. 
The leftmost column lists the possible dominant effects in /i, 



Eq. (111. The parameter t* {t* > t) depends on A(0) and H 




< /3 < 1 /3 < < 1 


Step interaction 
Line tension 


Aoe-'^''^' Ao{l + cBaAot)-^ 
Aoy^l-t/f Ao{l-t/f) 



Our predictions, based on Eq. (12 1 with ADL kinetics 



can be extended to TDL kinetics. The mobility M then 



reduces to Thus, we obtain (|13|) or (|14|, regardless 



of a, for step-interaction or line-tension dominated /i. 

Next, we compare our predictions with observations of 
Si(OOl) Oil] and Ag(llO) [8] corrugations. In Si(OOl), 
with i = 2DJk > 1000 nm [TD] and terrace width 



At 



< 



10 nm m El [TU], g(|V/i|) 



Aw 



> 100 which 



suggests ADL kinetics. We find decay laws comparing 
(i) the kinetic factor P, (3 ~ ^ < 0.01, with the aspect 
ratio squared, 

■ 'y 

of step interactions, ^, with (^)^. In 5 a ss 10^"^ 



(4^)^; and (ii) the relative strength 

Ay 

93 
91 

and thus /3 > a^. Also, « 1/30 and |^ > 1 |26], and 
thus 1^ ^ (^)^- Equation (14 1 follows, in agreement 
with the decay in [5]. In [7] a ^ 10~^, ly « 1/15 and 
w 100 So, /3 < a2 and ^ > )2. Equation (jis]) 



^) and 



follows, in agreement with the inverse linear decay in^7[. 

We now discuss observations of Ag(llO) [S] where 
step interactions are mainly entropic [27]. By 
f3 < 10-3 [8 and a « 1/15 [28,, we have l3 < 
We estimate by oi = ^ - ^ ln(coth 

93 « "'°°>^ [sinh(^)]-^ [2], where et is the kink 
formation energy, 0.04 eV< ek < 0.1 eV 0153], a = 1.4 
A, ao « 4 A, and T = 210 K; thus, < ^ < I. With 

ly « 2/25 [HI [3D], 1^ = O(^); thus, our criterion for 
step energetics appears inconclusive for scaling. Possible 
reasons are deviations of initial profiles from sinusoidal 
ones and anisotropies in Ag(llO), for which the model 
in [24j may be relevant. Although further study of the 
dynamics with reliable boundary conditions at facets is 
suggested, we view the condition /3 < as an indicator 
of evolution toward inverse linear decay [Hj. 

Our work forms a basis for a general approach to mor- 
phological evolution below roughening. Extensions in 2D 
include the ES barrier, long-range step interactions, step 
edge diffusion, anisotropy of step stiffness, and material 
deposition. Inclusion of the ES barrier [31] with rates 
and fed amounts effectively to k = 2{l/ku + l/kd)~^ in 
Eq. (12) [52]. Step-edge diffusion contributes to longitu- 



dinal fluxes but may not be important for Si(OOl), where 
ADL kinetics can dominate [2]. Anisotropic terrace dif- 



fusion, which is present in Si(OOl) and Ag(llO), is not 
expected to alter the main decay laws presented here. 
Connections of initial conditions and solutions for 



Eq. (12 1 to actual experimental situations have yet to 
be explored. Our scaling Ansatz should be tested for 
realistic initial profiles. Despite mode coupling [20] 



caused by the nonlinear PDE (12 1, our scaling should be 
valid for a range of prevailing wavelengths [T] [HI [20] . 

Other predictions of our approach include crossovers 
from exponential to inverse linear profile decay via 
aspect-ratio changes of the surface shape. Our work 
should stimulate further studies and relaxation exper- 
iments on surfaces below rougnening. 
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